---
title: "CARRS RDD analysis using a standardized cutoff - pair wise results"
knit: (function(input, ...) {rmarkdown::render(input, output_file="")})
output: pdf_document
geometry: "left=2cm,right=2cm,top=1cm,bottom=1cm"
date: "`r format(Sys.time(), '%d %B, %Y')`"
---

```{r , setup, include=FALSE, eval = TRUE}
knitr::opts_chunk$set(
   comment = '', warning=FALSE,  message=FALSE, results= FALSE , fig.width = 9,fig.height = 7, echo=FALSE, root.dir =""
)

 knitr::opts_knit$set(
   root.dir = ""    
 )
```

```{r packages}
rm(list=ls())
  library(tidyverse)
  library(gridExtra)
  library(dplyr)
  library(rdrobust)
  library(rddensity)
  library(viridisLite)
  library(foreign)
  library(ggplot2)
  library(data.table)
  library(doBy)
  library(socviz)
  library(openxlsx) 
  library(modelsummary)
  library(labelled)
  library(gtsummary)
  library(gt)
  library(kableExtra)
```

```{r loaddata}
data<-data.frame(read.csv("carrs_bp.csv"))
data<-as.data.table(data)
```

##calculating the bandwidth for each sub group

```{r}
# Males and females separately
rv.update = function(a,b,c) {

    out=c[a==1 & female==b,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                             (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
a_b<<-out

}

rv.update('w02',1,data)
w02_female<-a_b[w02==1 & female==1,]
rm(a_b)

rv.update('w02',0,data)
w02_male<-a_b[w02==1 & female==0,]
rm(a_b)

rv.update('w24',1,data)
w24_female<-a_b[w24==1 & female==1,]
rm(a_b)

rv.update('w24',0,data)
w24_male<-a_b[w24==1 & female==0,]
rm(a_b)


# total population
rv.updateMT = function(a,c) {
  out=c[a==1,runningvar:=ifelse((((sys-140)/(sd(sys))) == ((dias-90)/(sd(dias)))),
                                             (sys-140)/(sd(sys)),
                                            ifelse((((sys-140)/(sd(sys))) > ((dias-90)/(sd(dias)))),((sys-140)/(sd(sys))),((dias-90)/(sd(dias)))))]
a_b<<-out
}

rv.updateMT('w02',data)
w02_total<-a_b[w02==1 ,]
rm(a_b)

rv.updateMT('w24',data)
w24_total<-a_b[w24==1 ,]
rm(a_b)
```
##RDD Outputs function

```{r fuc}
  carrs_rdd = function(a,b,d) {
    d$Y<-a
    d$X<-b

    out =  rdrobust(y = d$Y, d$X, kernel = "triangular", p = 1, 
                    c = 0, rho = 1 , all = T)
    summary(out)
    
    ret <- data.frame(
sex = c(formatC(round(out$coef[1,1],2),2,format="f"), formatC(round(out$pv[3,1],2),2,format="f"), 
            paste0( "(", formatC(round(out$ci[3,1],2),2,format="f"), " - ", formatC(round(out$ci[3,2],2),2,format="f"), ")" ), 
            formatC(round(out$bws[1,1],2),2,format="f"), formatC(out$N_h[1] + out$N_h[2]))
  )
    rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
    
    output<<-ret %>% add_rownames  
  }
```

##Calling the specific subroups and saving the output

```{r}
rdd = function(a) {
  carrs_rdd(a$sysdiff,a$runningvar,a)
}

rdd(w02_male)
w02_male_sys <-output %>% mutate(Male=sex)%>%select(-sex)

rdd(w24_male)
w24_male_sys <-output%>%mutate(Male=sex)%>%select(-sex)

rdd(w02_female)
w02_female_sys <-output%>%mutate(Female=sex)%>%select(-sex)

rdd(w24_female)
w24_female_sys <-output%>%mutate(Female=sex)%>%select(-sex)

rdd(w02_total)
w02_total_sys <-output%>%mutate(Total=sex)%>%select(-sex)

rdd(w24_total)
w24_total_sys <-output%>%mutate(Total=sex)%>%select(-sex)
```

```{r diastolic}
rdd = function(a) {
  carrs_rdd(a$diasdiff,a$runningvar,a)
}

rdd(w02_male)
w02_male_dias <-output%>%mutate(Male=sex)%>%select(-sex)

rdd(w24_male)
w24_male_dias <-output%>%mutate(Male=sex)%>%select(-sex)

rdd(w02_female)
w02_female_dias <-output%>%mutate(Female=sex)%>%select(-sex)

rdd(w24_female)
w24_female_dias <-output%>%mutate(Female=sex)%>%select(-sex)

rdd(w02_total)
w02_total_dias <-output%>%mutate(Total=sex)%>%select(-sex)

rdd(w24_total)
w24_total_dias <-output%>%mutate(Total=sex)%>%select(-sex)
```

```{r}
table_results = function(x,y) {
  x<-x
  y<-y
  c<-merge(x,y,by="rowname")
  c<-c%>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))
out<<-c
}
```


```{r}
table_results(w02_total_sys,w24_total_sys) 
pairs_table_tdiag<-out
rm(c)
colnames(pairs_table_tdiag)[2] <- "Pair A" 
colnames(pairs_table_tdiag)[3] <- "Pair B"   
colnames(pairs_table_tdiag)[1] <- " " 

table_results(w02_female_sys,w24_female_sys)
pairs_table_fdiag<-out
rm(c)
colnames(pairs_table_fdiag)[2] <- "Pair A" 
colnames(pairs_table_fdiag)[3] <- "Pair B"   
colnames(pairs_table_fdiag)[1] <- " " 

table_results(w02_male_sys,w24_male_sys)
pairs_table_mdiag<-out
rm(c)
colnames(pairs_table_mdiag)[1] <- " " 
colnames(pairs_table_mdiag)[2] <- "Pair A" 
colnames(pairs_table_mdiag)[3] <- "Pair B"   

comb_sys_table<-rbind(pairs_table_tdiag,pairs_table_fdiag,pairs_table_mdiag)


table_results(w02_total_dias,w24_total_dias) 
pairs_table_ttreat<-out
rm(c)
colnames(pairs_table_ttreat)[2] <- "Pair A" 
colnames(pairs_table_ttreat)[3] <- "Pair B"   
colnames(pairs_table_ttreat)[1] <- " " 

table_results(w02_female_dias,w24_female_dias)
pairs_table_ftreat<-out
rm(c)
colnames(pairs_table_ftreat)[2] <- "Pair A" 
colnames(pairs_table_ftreat)[3] <- "Pair B"   
colnames(pairs_table_ftreat)[1] <- " " 

table_results(w02_male_dias,w24_male_dias)
pairs_table_mtreat<-out
rm(c)
colnames(pairs_table_mtreat)[1] <- " " 
colnames(pairs_table_mtreat)[2] <- "Pair A" 
colnames(pairs_table_mtreat)[3] <- "Pair B"   

comb_treat_table<-rbind(pairs_table_ttreat,pairs_table_ftreat,pairs_table_mtreat)
comb_pair_table<-cbind(comb_sys_table,comb_treat_table)
comb_pair_table<-comb_pair_table[ , -c(4) ]

colnames(comb_pair_table) <-c('','Pair A','Pair B',  'Pair A','Pair B')
```

#combining systolic and diastolic
```{r }
  filename <- paste0("pairwise_combined_bpdiff.tex")
  caption=paste0("Effect of home-based screening on systolic and diastolic blood pressure by survey pair")

x = kable(head(mtcars), "html")
combined<-kbl(comb_pair_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccccccc", format = "latex",digits = 2) %>%
    add_header_above(c(" "=1, "Systolic" = 2, "Diastolic" = 2)) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Female", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10) %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Male", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15) %>%
    footnote(general = c("Pair A: Wave 1 and 3 and Pair B: Wave 3 and 5. The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect for individuals with a standardized blood pressure of zero and CIs resulting from robust bias-corrected standard errors clustered at the individual-level.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")
```

